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ABSTRACT 


Without  nuclear  testing,  advanced  simulation  and  experimental  facilities,  such  as 
the  National  Ignition  Facility  (NIF),  are  essential  to  assuring  safety,  reliability,  and 
effectiveness  of  the  nuclear  force;  these  capabilities  are  invaluable  to  the  nation’s 
Stockpile  Stewardship  Program  (SSP).  A  significant  information  gap  exists  in  the 
hydrodynamic  response  to  nuclear  detonations  that  occur  near  the  earth’s 
surface.  Numerical  simulation  methods  were  used  to  evaluate  the  hydrodynamic 
response  of  earth-like  materials  and  to  develop  the  energy  coupling/partitioning 
curve  for  low  yield  nuclear  detonations  close  to  the  earth’s  surface.  Using  LLNL’s 
supercomputers  and  GEODYN  hydrodynamic  code,  the  properties  of  stress, 
pressure,  and  energy  were  evaluated  for  twelve  simulated  2.5kT  detonations;  six 
above  the  surface  and  six  below  the  surface.  The  results  indicate  stronger  air 
blasts  for  detonations  above  or  near  the  surface  and  that  energy  coupling  into  the 
ground  changes  rapidly  with  detonation  location  over  a  very  small  range  between 
the  above-ground  and  below-ground  interface.  This  work  serves  to  provide  a 
baseline  model  to  evaluate  stress,  pressure,  and  energy  in  relation  to  nuclear 
yield  close  to  the  earth’s  surface.  The  results  support  a  better  understanding  of 
the  physics  of  near-surface  detonations  and  also  assist  in  planning  future 
experimental  work  at  NIF. 
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I.  INTRODUCTION 


From  July  1945  to  September  1992,  the  United  States  conducted  1,054 
nuclear  tests  and  executed  two  nuclear  attacks.  Considering  the  tests  performed 
by  other  nuclear  weapon  states,  the  number  of  nuclear  tests  done  by  the  U.S. 
nearly  doubles  that  number.  Since  the  initial  scientific  discovery  of  nuclear 
fission  and  the  early  development  of  nuclear  weapon  technology,  the 
international  situation  has  matured  considerably,  and  there  has  been  an 
increasing  concern  about  the  potentially  harmful  environmental  and  public  health 
effects  of  nuclear  testing,  in  addition  to  the  devastating  results  of  nuclear  war. 
Over  the  decades,  there  have  been  many  efforts  on  the  part  of  the  United  States 
to  reduce  the  threat  of  nuclear  weapons.  These  efforts  culminated  in  August  of 

1995  when,  then  President  Bill  Clinton,  announced  the  Comprehensive  Test  Ban 
Treaty  (CTBT),  which  would  institute  a  worldwide  ban  on  the  testing  of  nuclear 
weapons  of  any  yield.  This  treaty  was  signed  by  the  United  States  in  September 

1996  at  the  United  Nations,  though  the  agreement  has  not  yet  been  ratified  by 
the  U.S.  Senate  and,  thus,  has  not  officially  entered  into  force.  Nevertheless,  the 
U.S.  and  other  nations  have  complied  with  the  major  features  of  the  treaty,  and  a 
de  facto  ban  on  nuclear  weapons  testing  has  emerged. 

Given  the  current  de  facto  ban  on  nuclear  testing,  in  order  to  remain 
confident  in  the  safety,  reliability,  security  and  deterrent  effectiveness  of  the 
existing  nuclear  weapons  arsenal,  while  retaining  proficiency  in  nuclear 
technology,  methods  other  than  nuclear  testing  have  been  developed  [1].  The 
Science-based  Stockpile  Stewardship  (SBSS)  Program  was  established  to: 
increase  the  understanding  of  the  current  nuclear  stockpile,  predict  and  evaluate 
potential  problems  as  these  nuclear  weapons  age,  fix  or  re-manufacture  weapon 
components,  and  maintain  proficiency  and  credibility  to  support  the  deterrent 
value  of  potential  nuclear  weapon  use  in  the  future  [2].  According  to  the  U.S. 
DOE/NNSA,  “Computer  models  and  advanced  experimental  capabilities  that 
provide  accurate  predictive  capability  of  weapon  performance  in  the  absence  of 
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nuclear  testing  became  the  main  goal  of  the  National  Nuclear  Security 
Administration’s  Stockpile  Stewardship  Program. ”[2]  The  National  Ignition  Facility 
(NIF)  at  Lawrence  Livermore  National  Laboratory  serves  as  a  vital  tool  to  the 
nation’s  Stockpile  Stewardship  Program.  NIF  is  currently  the  largest  laser 
system  ever  built  in  the  world.  It  has  the  capability  to  achieve  pressure  and 
temperature  conditions  that  are  replicated  only  in  nuclear  weapons,  the  sun,  and 
stars,  albeit  on  a  much  smaller  scale  [1]. 

During  the  era  of  nuclear  weapons  testing,  there  were  many  experiments 
conducted  that  measured  the  hydrodynamic  properties  and  shock  response  of 
materials  above  and  below  ground.  There  is  a  wealth  of  knowledge  documented 
from  the  nuclear  tests  that  were  conducted,  but  the  majority  of  these  nuclear 
tests  took  place  either  at  altitudes  well  above  ground  level  or  at  depths  of  burial 
well  below  the  ground  surface.  Consequently,  there  is  little  data  of  the 
hydrodynamic  response  of  materials  to  a  nuclear  detonation  at  or  near  the 
earth’s  surface.  Given  the  current  ban  on  nuclear  testing,  large-scale  testing  of 
weapons  cannot  fill  such  data  gaps;  however,  LLNL’s  National  Ignition  Facility 
offers  a  unique  capability  to  generate  data  in  this  region  of  interest  that  cannot  be 
duplicated  by  any  other  conventional  means.  Employing  the  laser  technology  of 
NIF,  along  with  advanced  numerical  computation  methods,  data  can  be  produced 
that  give  insight  to  both  the  hydrodynamic  response  of  materials,  and  the  shock 
physics  that  takes  place  during  and  after  a  nuclear  detonation  at  or  near  the 
earth’s  surface.  The  data  generated  will  serve  particular  interests  in  the  fields  of 
nuclear  detonation  detection  and  verification,  nuclear  forensics,  and  structural 
survivability.  The  primary  objective  of  this  thesis  research  is  to  perform  a  series 
of  simulation  calculations  to  analyze  and  support  the  planning  of  NIF  experiments 
designed  to  characterize  the  shock  response  of  a  prototypical  material  from  near¬ 
surface  detonations,  and  to  generate  an  energy  coupling  plot  that  shows  the 
air/ground  energy  coupling  and  the  energy  partitioning  between  air  and  ground 
that  results  from  nuclear  detonations  occurring  at  or  near  the  earth’s  surface. 
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II.  PROBLEM  SETUP 


A.  DEFINING  THE  PROBLEM 

The  principal  factors  that  distinguish  a  nuclear  detonation  from 
conventional  detonations  are  specific  energy  and  yield.  The  two  nuclear 
weapons  used  in  the  attacks  conducted  by  the  United  States  on  Hiroshima  and 
Nagasaki  had  yields  of  12-15  kilotons  of  TNT  equivalent  and  20-22  kilotons  of 
TNT  equivalent,  respectively.  In  a  NIF  experiment  simulating  an  actual  weapon 
detonation,  the  experimental  setup  must  be  planned  to  properly  scale  down  the 
physical  dimensions  and  characteristics  of  the  event.  In  order  to  properly  mimic 
the  conditions  of  a  low  yield  nuclear  explosion,  the  parameters  of  energy,  length, 
and  time  must  be  scaled  down  to  the  lab  frame.  The  National  Ignition  Facility 
has  192  laser  beams  each  with  the  capability  of  delivering  approximately  9.4kJ  of 
energy  for  a  combined  total  of  1.8MJ  of  light  energy.  NIF  can  generate  a  peak 
power  of  500  trillion  watts,  which  is  1,000  times  the  electric  generating  power  of 
the  entire  country  [3].  High-precision  lasers,  such  as  these,  allow  for  the 
achievement  of  specific  energies  that  are  near  those  of  nuclear  devices.  This 
opens  up  a  window  of  opportunity  to  investigate  the  physics  of  near-surface 
nuclear  detonations.  For  the  experiments  being  considered,  only  four  of  the  192 
laser  beams  will  be  used,  with  each  contributing  only  2.5kJ  (i.e.,  operating  at 
roughly  25%  capacity),  to  achieve  a  total  deposited  energy  of  lOkJ  focused  into  a 
2mg  cylindrical  hohlraum,^  which  yields  a  specific  energy  of  5MJ/g.  This 
cylindrical  hohlraum  has  a  height  and  diameter  of  about  2mm  and  a  volume  100 
times  smaller  than  that  of  a  pencil  eraser. 

B.  SCALING 

To  scale  the  experiment  down  to  the  laboratory  frame,  the  Glasstone/ 
Dolan  cube  root  scaling  law  was  used: 

A  small  gold  cylinder  named  after  the  German  word  for  “hollow  room”  that  serves  as  the 
target  for  a  NIF  experiment. 
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Theoretically,  a  given  pressure  will  occur  at  a  distance  from  an 
explosion  that  is  proportional  to  the  cube  root  of  the  energy  yield. 
Full-scale  tests  have  shown  this  relationship  between  distance  and 
energy  yield  to  hold  for  yields  up  to  (and  including)  the  megaton 
range.  [4] 

For  the  anticipated  experiment  relating  a  2.5kT  full-scale  test  to  a  NIF 
energy  deposition  of  10kJ,  a  scaling  factor  of  1000  can  be  derived  as  follows: 

{2.5kTI\mY^  =  [2.5^r* (4.1 84x1 0’ / 1 ~  1000 
where  lytr  =  4.1 84  *10‘V 

To  elaborate,  in  an  actual  2.5  kT  yield  detonation,  the  energy  release  is: 
2.5kT  =  10.46  X  10^^  J  =  10.46  x  10^  kJ 

This  experimental  NIF  energy  deposition  of  lOkJ  in  comparison  with  the 
detonation  release  of  2.5kT  represents  an  energy  reduction  factor  of: 

(2.5kT/10kJ)  =  (10.46  x  10^  kJ/10  kJ)  ~  10® 

Applying  Glasstone  and  Dolan’s  cube  root  scaling  factor  rule,  the  resulting 
scaling  factor  is  10®,  or  a  factor  of  1000  [4].  Another  way  to  characterize  the 
cube  root  scaling  factor  rule  as  applied  to  distance  scaling  is  to  recognize  that: 

Da=  Dr  *  (Wa/  Wr)"'® 

where: 

Da=actual  distance 
Dr=reference  distance 
Wa=actual  yield 
Wr=reference  yield 

Therefore,  in  the  case  being  considered  where  the  energy  reduction  factor 
is  10®,  the  distance  scaling  is  a  factor  of  one  thousand.  A  similar  calculation  is 
also  applied  to  the  time  parameter.  As  a  result,  with  both  time  and  distance 
properly  scaled,  the  physical  properties  such  as  particle  velocity,  pressure, 
density,  and  temperature  are  accurately  reflected  without  scale  factors.  The 
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primary  properties  of  interest  in  this  thesis  are  those  of  pressure  and  energy, 
which  will  be  used  to  construct  the  energy-coupling  curve  as  the  height  of  burst/ 
depth  of  burial  is  varied  for  a  2.5kT  simulated  explosion. 

For  example,  Figure  1  presents  the  results  of  a  simulation  of  a  10kJ  NIF 
experimental  explosion  at  a  height  of  burst  of  3cm  (corresponding  to  a  2.5kT 
nuclear  detonation  at  a  height  of  burst  of  30m),  in  which  the  pressure  wave 
reaching  a  maximum  of  .8MPa  (8  bar),  arrives  in  approximately  80  microseconds 
at  a  sensor  location  (Gauge  #2)  15  cm  from  the  explosion.  This  would  scale 
back  to  the  real  world  reference  frame  as  a  shock  with  the  same  peak  pressure 
of  8  bar  arriving  in  approximately  80  milliseconds  at  a  distance  of  150  meters 
away  from  an  explosion  of  2.5kT  yield  at  a  height  of  burst  of  30  meters. 
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III.  PHYSICAL  CONFIGURATION 

A.  NIF  CONFIGURATION 

For  a  series  of  anticipated  NIF  experiments  to  investigate  near-surface 
detonation  coupling,  four  of  the  192  laser  beams  will  be  used.  Each  beam 
contributes  2.5  kJ  of  energy  over  two  nanoseconds  that  will  be  deposited  in  a 
2mg  cylindrical  hohlraum.  The  four  lasers  will  combine  inside  the  near-vacuum 
target  chamber,  as  shown  in  Figures  2  and  3,  and  focus  directly  into  the  gold 
hohlraum  target.  The  gold  will  essentially  vaporize  and  expand  as  a  gold  gas  pill. 
This  is  similar  to  physical  conditions  in  the  detonation  of  a  full  scale  nuclear 
device  where  all  components  are  vaporized  and  contribute  to  the  ensuing  fireball. 
The  gold  hohlraum  will  be  located  30mm  above  the  surface  of  a  cylindrical  core 
of  a  ceramic  material  called  Macor2,  which  for  the  purpose  of  these  experiments, 
simulates  the  solid  material  of  the  earth’s  surface.  The  properties  of  the  Macor 
ceramic,  some  of  which  are  listed  Table  1,  are  similar  to  the  properties  of  many 
solid  earth  materials,  in  particular  limestone.  The  Macor  ceramic  in  the 
experiment  serves  as  a  model  for  the  ground  response  to  a  nuclear  detonation 
over  limestone. 


Density 

2520  kg/m^ 

Young’s  Modulus 

64.1  GPa 

Shear  Modulus 

25.4  GPa 

Poisson’s  Ratio 

0.26 

Quasi-static  compressive  strength 

345  MPa 

Table  1 .  Physical  Properties  of  Macor  glass  ceramic  (From  W.Chen,  1997, 

p.1310) 


2  Macor  is  a  machinable  white  glass  ceramic. 
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Figure  3.  Target  Chamber  Exterior 

The  experimental  target  consists  of  a  hohlraum  target  (representing  the 
scaled  down  nuclear  device)  and  the  pre-fabricated  Macor  cylinder  which  will 
include  three  pressure  sensors  embedded  into  the  structure  in  addition  to  the 
hohlraum,  and  two  additional  pressure  sensors  mounted  outside  of  the  Macor 
object  in  order  to  measure  pressure  in  the  region  representing  the  near-surface 
atmosphere.  From  the  focal  point  of  the  laser  beam,  the  pressure  gauges^  will 
be  located  at  the  straight-line  distances  shown  in  Table  2. 


3  PG1  and  PG2  represent  the  two  atmospheric  pressure  sensors  external  to  the  Macor  object 
while  GPG1  through  3  represent  the  pressure  sensors  embedded  within  the  Macor  object. 
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Pressure  Gauge  1  (PG1) 

10cm 

PG2 

15cm 

Ground  Pressure  Gauge  1  (GPG1) 

5cm 

GPG2 

7cm 

GPG3 

10cm 

Table  2.  Distances  from  target  to  pressure  sensors 


Pressure  Gauges  1  and  2 

r  '' 

PG1  PG2 
Ground  Pressure  Gauges  1. 2.  af>d  3 


GPG1 

f 

i 

Figure  4.  Set  Up  Drawing  After  [6] 

Figure  4  shows  a  diagram  of  the  experiment,  self  contained  in  a  cylindrical 
casing  that  will  be  inserted  on  the  diagnostic  instrument  manipulator  (DIM) 
depicted  in  Figure  5.  The  subsurface  sensor,  GPG2,  has  a  quartz  piezoelectric 
element  with  a  dynamic  range  of  0-2  GPa  and  a  bandwidth  of  1-100  MFIz.  The 
subsurface  sensors  GPG1  and  GPG3,  located  at  45-degree  angles  from  the 
laser  vector  in  its  focal  point,  have  polyvinylidine  fluoride  (PVDF)  piezoelectric 
elements,  which  provide  longer  recording  times,  and  each  has  a  dynamic  range 
of  0-2  GPa  and  a  bandwidth  of  1-100  MFIz.  PVDF  is  a  piezoelectric  plastic 
polymer  that  is  commonly  used  as  a  shock  sensor  [7],  The  above-surface 
sensors  PG1  and  PG2  have  PVDF  elements  and  have  a  dynamic  range  of  0-6 
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4  converging 
NIF  beams 


MPa  and  a  bandwidth  of  30kHz  to  3MHz.  The  Macor  cylinder  will  be  inserted  into 
the  NIF  target  chamber  attached  to  the  diagnostic  instrument  manipulator  (DIM), 
shown  in  Figures  5  and  6. 


Figure  5.  DIM  and  Chamber 


Figure  6.  Close  Up  View  of  the  DIM 


B.  CONFIGURATION  OF  NIF  EXPERIMENTAL  DESIGN  AND 
CORRESPONDING  GEODYN  SIMULATIONS 

In  order  to  properly  model  the  experiment  simulating  a  nuclear  detonation 
over  or  in  limestone,  a  combination  of  materials  with  similar  properties  to  that  of 
limestone  and  air  must  be  used.  With  Macor  as  the  solid  material,  the  gas  mix 
was  selected  by  empirical  methods  to  mimic  as  close  as  possible  atmospheric 
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density.  The  required  gas  mix  was  calculated  to  consist  of  65%  Ne,  20%  Ar,  10% 
Kr,  and  5%  Xe.  This  mixture  yields  similar  results  in  terms  of  transmission  and 
absorption  of  the  soft  x-rays  emitted  in  order  to  produce  a  scaled  deposition  of 
the  x-rays  and  the  “double  peak”  intensity  plot  that  is  characteristic  of  nuclear 
detonations.  The  first  peak  in  light  intensity  is  due  to  the  initial  reaction  and 
expansion  of  hot  gas  in  the  form  of  a  fireball.  Quickly,  the  shock  wave  from  the 
explosion  expands  past  the  fireball  and  the  air  outside  the  fireball  is  compressed 
by  the  shock  wave  to  the  point  of  becoming  opaque  to  light.  As  the  shock  wave 
expands  spherically,  the  pressure  front  decreases  and  the  light  from  the  fireball  is 
once  again  visible  [8].  Table  3  presents  a  summary  comparison  of  the  properties 
of  the  Macor/artificial  atmosphere,  and  limestone/normal  atmosphere  which 
demonstrates  their  similarities. 


PROPERTY 

GAS  MIX  OVER  MACOR 

AIR  OVER  LIMESTONE 

‘Atmospheric’  density 

.00153  g/cc 

.001204  g/cc 

‘Ground’  density 

2.52  g/cc 

2.59  g/cc 

‘Ground’  sound  speed 

4.58  km/s 

3.67  km/s 

Table  3.  MACOR  and  Limestone  Comparison  After  [6] 


To  achieve  the  objective  of  this  thesis,  a  series  of  simulations  using 
LLNL’s  GEODYN4  computer  code  was  carried  out  to  analyze  the  experimental 
configuration  for  planning  purposes  and  also  to  develop  response  projections  to 
establish  the  ground-coupling  curve.  The  first  set  of  simulations  were  carried  out 
to  represent  the  hydrodynamic  response  to  detonations  (and  NIF  experiments 
representing)  above-ground  detonations.  The  second  set  treated  below-ground 
detonations. 

1.  Above  Ground  Configuration 

Figure  7  presents  an  experimental  plan  and  the  corresponding  simulation 
layout  for  the  above-ground  case. 

^  More  information  about  GEODYN  can  be  found  at  https://e-reports- 
ext.llnl.gov/pdf/326532.pdf. 
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Figure  7.  3cm  Height  of  Burst  Calculation  Set  Up;  Z-and  R-  axis  in 

millimeters 

In  this  portion  of  the  simulations,  a  simplified  2D  axisymmetric  geometrical 
representation  of  the  physical  configuration  was  constructed.  Figure  7  shows  the 
output  of  LLNL’s  graphics  software,  VisIT^,  at  time  zero,  as  the  initial  set  up  in  the 
GEODYN  computer  modeling  code.  In  the  above-ground  model  setup,  there  are 
five  materials — Macor,  the  atmospheric  gas-mix,  gold,  “air,”®  and  steel.  The 
purposes  of  the  first  three  materials  are  obvious,  but  the  latter  two  require 


5  VisIT  is  an  interactive  parallel  visualization  and  graphical  analysis  tool  for  scientific  data 
capable  of  handling  data  on  the  terabyte  scale. 

®  The  “air”  is  given  an  extremely  low  density  to  simulate  the  vacuum  conditions,  through 
which  the  lasers  propagate,  inside  the  NIF  target  chamber. 
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additional  explanation.  Within  the  computer  model,  there  must  be  a  replication  of 
the  vacuum  barrier  through  which  the  lasers  propagate  to  deposit  the  energy 
inside  of  the  gas  mix  for  an  above-ground  calculation.  In  this  simulation,  a  hollow 
steel  cone  is  used  to  represent  the  vacuum  barrier  and  to  keep  that  vacuum 
separate  from  the  gas  mix.  (See  Appendix  A)  Essentially,  the  simulation  begins 
after  all  the  energy  from  the  laser  is  deposited  into  the  gold  hohlraum.  This  was 
achieved  by  creating  a  steel  cone  with  an  apex  at  the  laser  focal  point.  In  the 
simulation,  the  steel  cone  is  filled  with  “air,”  but  given  a  density  of  almost  zero,  so 
that  it  does  not  contribute  to  the  forming  shock  wave.  The  actual  pressure 
anticipated  inside  of  the  NIF  target  chamber  is  close  to  10'^  millitorr  at  the  time  of 
a  target  shot.  By  changing  the  location  of  the  apex  of  the  steel  cone  (and  the 
attached  hohlraum  target),  the  desired  height  of  burst  for  the  calculation  is 
achieved.  In  order  to  record  the  changes  in  the  hydrodynamic  properties  as  the 
experiment  evolves,  marker  files  were  embedded  throughout  the  configuration 
with  the  primary  markers  being  geometrically  located  at  the  same  places  that  the 
subsurface  and  aboveground  pressure  sensors  are  located.  Each  of  the  five 
materials  is  indexed  for  identification  with  label  values  from  0  to  4  and  given  a 
density  and  a  strength  parameter  and  defined  by  a  series  of  coordinate  points. 
The  input  file  for  the  10mm  height  of  burst  simulation  can  be  found  in  Appendix  A 
for  further  clarification. 

2.  Below  Ground  GEODYN  Configuration 

The  depth  of  burial  simulation  set  up  is  similar  to  that  of  the  height  of  burst 
simulations,  with  the  exception  that  there  is  no  need  to  incorporate  the  steel 
cone,  as  the  energy  is  deposited  directly  into  a  hohlraum  embedded  in  the  Macor 
block.  Figure  8  shows  the  initial  set  up  of  the  simulation,  in  VisIT,  for  the  below¬ 
ground  case.  The  subsurface  calculations  are  done  from  an  idealized  point  of 
view.  The  Macor  is  uniformly  manufactured  and  the  Macor  representation  is 
uniform.  This  is  a  limitation  in  terms  of  representing  real-world  earthen 
structures,  which  are  far  from  uniform.  The  possible  addition  of  cracks,  porosity, 

or  other  heterogeneities  within  the  simulated  Macor  is  not  addressed  in  this  work 
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but  may  be  a  topic  for  future  study.  Additionally,  in  terms  of  mimicking  a  setup 
for  a  NIF  experiment  there  must  be  a  conduit  through  which  the  energy,  from  the 
laser,  is  deposited.  Essentially  this  would  involve  constructing  an  inverted  hollow 
cone  with  an  apex  beneath  the  surface.  The  amount  of  energy  that  remains  in 
the  region  of  the  cone  as  opposed  to  being  contained  completely  in  the  “ground” 
must  also  be  taken  into  consideration  when  calculating  the  energy  coupling.  It 
must  be  stressed  that  the  depth  of  burial  simulations  conducted  in  this  work  were 
done  from  an  ideal  standpoint.  This  would  be  the  real  world  equivalent  of  digging 
a  cavity  into  the  ground,  placing  a  10kJ  nuclear  device  in  the  cavity,  and  then 
covering  the  cavity  and  packing  it  to  the  ground’s  original  density.  These  depths 
of  burial  calculations  serve  as  a  baseline  used  to  construct  the  energy  coupling 
curve. 
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IV.  PHYSICS  BEHIND  THE  CODE 


A.  GEODYN  BASICS 

In  order  for  a  hydrodynamic  simulation  to  have  validity,  it  must  be 
constrained  so  as  to  obey  the  laws  of  physics.  The  hydrodynamic  code  used  in 
this  thesis,  called  GEODYN,  was  developed  at  LLNL  and  incorporates  physical 
models  to  fully  describe  a  broad  range  of  phenomena  including  shock  and 
thermodynamic  behavior.  It  is  an  Eulerian  code  with  adaptive  mesh  refinement 
(AMR).  There  are  two  basic  types  of  dynamic  codes,  Lagrangian  and  Eulerian. 
A  Lagrangian  code’s  mesh,  or  background,  moves  with  the  material,  so  no  mass 
flows  between  cells.  In  an  Eulerian  code,  the  mesh  or  background  is  stationary 
and  the  material  is  allowed  to  move  though  stationary  cells.  An  Eulerian  code  is 
analogous  to  looking  at  a  dust  storm  through  a  window  whereas  a  Lagrangian 
code  is  analogous  to  floating  along  with  an  individual  dust  particle  in  the  storm. 
The  adaptive  mesh  means  that  the  code  has  the  ability  to  vary  the  level  of  detail 
of  the  background.  In  keeping  with  the  dust  storm  analogy,  the  adaptive  mesh  is 
analogous  to  having  the  ability  to  simultaneously  have  small  microscopic 
windows  to  capture  the  details  of  small  grains  and  larger  windows  for  use  when 
detail  is  not  needed.  An  Eulerian  code  with  adaptive  mesh  refinement  such  as 
GEODYN  allows  for  rigorous  high  numerical  resolution  in  areas  in  one  part  of  a 
problem  and  less  refined,  less  sensitive  areas  in  another. 

When  dealing  with  shock  physics,  high  temperatures,  and  nuclear 
interactions,  the  laws  of  physics  that  constrain  the  results  are  numerous.  The 
two  primary  areas  of  physics  that  are  of  concern  are  thermodynamics  and 
continuum  mechanics.  The  GEODYN  code  simplifies  the  matter  without 
deviation  from  real  experimental  results  by  implementing  the  various  physics 
constraints  into  a  formulation  based  on  both  the  1st  and  2nd  laws  of 
thermodynamics,  and  the  laws  of  conservation  of  mass,  momentum  and  energy 
[9]. 
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GEODYN  is  a  complex  modeling  code.  In  order  to  provide  some 
background  information  on  key  physics  processes  taking  place  in  GEODYN,  the 
remaining  sections  in  this  chapter,  (Sections  B  and  C),  provide  some 
mathematical  derivations  of  some  of  the  physics  involved  in  terms  of 
thermodynamics  and  continuum  mechanics  respectively. 

B.  THERMODYNAMICS 

A  major  portion  of  GEODYN,  summarized  in  Table  4,  are  the  constitutive 
equations  that  express  specific  entropy,  specific  internal  energy,  stress,  and  heat 
flux  in  terms  of:  temperature  (0),  change  in  temperature  (A0),  deformation 
gradient  (F)  and  the  change  in  deformation  gradient  (dF). 


Specific  entropy  (s) 

s=s(0,A0,F,dF) 

Specific  internal  energy  (u) 

u=u(0,A0,F,dF) 

Stress  (T) 

T=T(0,A0,F,dF) 

Heat  flux  (q) 

q=q(0,A0,F,dF) 

Table  4.  Constitutive  Equations 


These  constitutive  equations  can  be  derived  from  the  Clausius-Duhem 
inequality^,  which  has  been  derived  below  for  clarification. 

Consider  the  following  five  parameters  of  interest: 

1 .  Specific  entropy  s 

2.  Mass  density  of  the  body  p 

3.  Internal  heat  supply/  (mass-time)  r 

4.  Outward  heat  flux  vector  q 

5.  Temperature  0 

^  The  Clausius-Duhem  inequality  is  used  in  continuum  mechanics  to  express  the  second  law 
of  thermodynamics. 
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The  entropy  (due  to  heat  transfer  into  a  system)  is  given  by  the  following 
formula^: 


ds 

dt 


=  j  pr  !  Odv  + 

Volume 


j  -q  !  Oda 

Surface 


where  dv  and  da  represent  the  defferentials  of  volume  and  area  respectively. 


However,  when  dealing  with  real,  irreversible  processes  the 
thermodynamic  equations  become  inequalities.  So,  the  change,  with  respect  to 
time  of  the  entropy,  would  be 

—  j  spdv>  j  rpl9dv+  j  -qlOda. 

Volume  Volume  Surface 


Using  the  divergence  theorem  (integral  over  area  is  equal  to  the  integral  of 
the  divergence  over  a  volume),  the  expression  becomes: 

—  j  spdv>  j  rplOdv-  j  V»(q/d)dv. 

Volume  Volume  Volume 


This  is  the  integral  form  of  the  Clausius-Duhem  inequality.  The  Clausius- 
Duhem  inequality  in  its  most  recognized  form  is  in  the  differential  form: 


—  >rlO-\lp^y*{qlO). 
dt 


This  formula  states  that  the  rate  of  entropy  change,  with  respect  to  time,  in  a 
system  is  greater  than  the  rate  of  entropy  change  due  to  heat  transfer  into  the 
system. 


C.  CONTINUUM  MECHANICS 

The  next  major  portion  of  GEODYN  deals  with  how  a  material  deforms,  or 
flows  over  time,  due  to  energy  input;  all  while  conserving  mass,  energy,  and 


®  The  mathematical  derivations  that  follow  are  summaries  based  upon  detailed  presentations 
in  [9]. 
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momentum.  The  equations  used  are  often  in  the  form  of  the  convective 
derivative  (deformation/  flow  equation);  for  example 


Ut  at 


This  equation  states  that  the  total  change,  with  respect  to  time,  is  equal  to 
the  change  in  the  fixed  frame,  plus  the  change  in  the  moving  frame.  In  effect, 
this  is  a  continuity  equation.  The  parameter  denoted  by  the  #  sign  represents  a 
history-dependant  scalar  parameter  such  as  porosity,  temperature,  or  plastic 
strain  and  V  (capitalized)  represents  the  velocity  of  the  fluid  medium.  The 
derivation  of  the  generic  continuity  equation  is  shown  below  for  the  purpose  of 
providing  clarification  for  what  is  taking  place  mathematically  in  the  GEODYN 
code. 


Let  L  be  some  intensive  (mass  independent)  property  defined  over  a 
control  volume  Q.  The  parameter  Q  represents  sources  and  sinks  in  the  fluid. 
The  rate  of  change  of  L  is  given  by  the  formula  below: 


d_ 

dt 


j  Ldv 

Q. 


-  (^  L  Vnda  -  j  Qdv 

Surface  Q 


where  dv  and  da  represent  the  defferentials  of  volume  and  area  respectively 


Applying  the  divergence  theorem  to  the  1st  term  on  the  right  side  of  the 
equal  sign  gives: 


A 

dt 


J  Ldv  =  -J  V  •  {LV)dv  -  J  Qdv 
0.0. 


Combining  the  terms  on  side  of  the  equation  gives: 

J(^  +  V.(TF)  +  0Jv  =  O 

The  above  condition  must  be  true  for  all  dv;  therefore,  the  following 
generic  continuity  equation  results: 
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+  V.(LF)  +  0  =  O 

at 


This  equation  states  that  the  amount  of  a  conserved  quantity,  inside  a 
defined  region,  can  only  change  by  the  amount  that  passes  through  the  boundary 
of  the  region  in  addition  to  any  internal  sources  or  sinks.  The  conservation  of 
mass,  energy,  and  momentum  all  follow  from  the  above  generic  continuity 
equation. 

In  GEODYN,  solid  is  modeled  as  an  extension  of  a  fluid,  with  the  added 
appropriate  strength  parameters.  A  LLNL  Mie-Gruneisen  equation  of  stated  was 
used  to  calculate  the  pressure  in  the  Macor  for  this  simulation.  For  the  gold 
source,  a  tabular  equation  of  state  provided  by  LLNL  was  used. 


^  The  Mie-Gruneisen  equation  of  state  is  often  used  to  measure  the  pressures  in  shock 
compressed  solids. 
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V.  METHOD 


A  series  of  simulation  calculations  were  performed  in  which  the  height  of 
burst  and  depth  of  burial  ranged  from  5mm  to  50mm.  A  simulation  was  run  at 
5mm,  10mm,  20mm,  30mm,  40mm,  and  50mm  for  both  above-ground  and 
below-ground  cases  for  a  total  of  12  simulations.  On  average,  each  simulation 
required  a  total  24  hours  of  computation,  broken  into  three  8-hour  periods  on 
LLNL  supercomputers.  Each  computation  was  considered  complete  once  the 
shock  front  passed  the  location  of  the  most  distant  shock  sensor,  which  on 
average  tended  to  be  about  100  microseconds  (real  time),  or  about  24  hours  of 
computation  time.  To  achieve  12  complete  data  sets  ranging  from  a  depth  of 
burial  of  5cm  to  a  height  of  burst  of  5cm  required  a  time  period  of  three  months. 
After  each  successful  calculation,  the  data  were  analyzed  using  both  MATLAB'io 
and  LLNL’s  visualization  software  VisIT'iT  GEODYN  has  the  ability  to  track 
various  parameters  during  the  experiment.  The  parameters  of  interest  are: 
pressure,  density,  material  type,  particle  velocity  in  the  r  and  z  coordinate 
directions,  sound  speed,  internal  energy,  total  energy,  and  temperature.  As 
stated  earlier,  in  order  to  construct  data  showing  the  energy  coupling  between 
ground  burst  and  air  shock,  the  primary  parameters  of  interest  are  pressure, 
particle  velocity,  time  and  total  energy.  Figures  9  through  13  present  the 
pressure  wave  results  for  the  10mm  height  of  burst  calculation  at  the  five  sensor 
positions. 


MATLAB  is  a  scientific  software  application  used  for  matrix  manipulation  and  plotting  that 
can  interface  with  other  programming  languages  such  as  C,  C++,  and  FORTRAN.  See 
http://www.mathworks.com/products/matlab/. 

"I  "I  FOR  ADDITIONAL  DESCRIPTION  OF  VISIT,  See 
https://wci.llnl.gov/codes/visit/about.html. 
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Height  of  Burst  Plots 


Figure  9.  Pressure  Wave  at  Air  Gauge  1  for  1cm  HOB 


Figure  10.  Pressure  Wave  at  Air  Gauge  2  for  1cm  HOB 
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Pr*siu>«  at  Ground  Pr*a$ui*  Gtug*  1  1cm  HOB 


Figure  1 1 .  Pressure  Wave  at  Ground  Gauge  1  for  1  cm  HOB 


Pressure  at  Ground  Pressure  Gauge  2  1cm  HOB 
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Figure  12.  Pressure  Wave  at  Ground  Gauge  2  for  1cm  HOB 
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Pressure  at  Ground  Pressure  Gauge  3  tcm  HOB 


Figure  13.  Pressure  Wave  at  Ground  Gauge  3  for  1cm  HOB 

The  most  useful  information  gathered  from  these  plots  is  the  peak 
pressure  and  the  time  of  arrival  of  the  peak  pressure  wave.  This  information 
served  as  a  check  for  the  calibration  of  the  PVDF  shock  sensors.  For  Figures  10 
through  13,  after  the  initial  pressure  waves  any  other  peaks  are  the  results  of 
boundary  reflections.  In  the  case  of  Figure  10  for  PG2  a  second  pressure  peak 
is  visible  at  120  microseconds.  This  second  peak  is  due  to  the  reflection  of  the 
initial  blast  wave  from  the  Macor/gas-mix  boundary.  For  Figures  11  though  13 
the  many  reflection  stress  waves  are  the  result  of  free  surface  reflections  from 
the  Macor/gas-mix  boundary.  Some  of  the  interesting  phenomena  shown  in  the 
plots  above  were  made  clear  when  the  data  were  analyzed  using  LLNL’s 
visualization  software,  VisIT.  For  example,  in  the  ground  sensors,  the  stress 
waves  in  tension  were  easily  identified  by  stepping  through  the  calculation  frame- 
by-frame  in  time.  In  the  simulation,  the  Macor  is  surrounded  by  a  free  surface 
and  thus  tension  waves  occur  at  the  intersection  between  the  Macor  block  and 
the  gas  mix.  Figure  14  shows  a  movie  frame  of  the  detonation  at  11 
microseconds,  with  pressure  on  the  left  and  material  on  the  right,  where  the  initial 
stress  wave  propagating  through  the  Macor  object  is  visible.  Figure  15  shows  a 
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movie  frame,  at  24  microseconds,  where  the  stress  waves  have  reached  the 
edges  of  the  Macor  object  and  reflected  tension  waves  are  clearly  visible. 


DB:  Header 

Cvcie;  -2147403047  Tlme:l  1  £084 

Pseudocolor 

Varcobi 

4.0X1 


I 


^"-O.OCD 
Max:  4.000 
Mh:  0.000 


Pseudocolof 


-oxitcieo 


^^■--a022S7 
Mate  0.1349 
Min: -0.02257 


-50 


-100 


Figure  14.  Stress  wave  propagating  through  the  Macor  object 
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Figure  15.  Tension  waves  propagating  through  the  Macor  object 
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Since,  in  the  actual  experiment,  the  Macor  test  object  will  be  held  on  the 
DIM,  the  real  fluctuation  in  pressure  can  be  expected  to  be  substantially  less. 
The  recovery  period  of  the  pressure  sensor  will  also  eliminate  much  of  the 
reflected  response,  as  the  sensors  will  not  be  able  to  record  continuously  as  in 
the  simulation. 

Depth  of  Burial  Plots 

The  pressure  waves  for  the  five  shock  sensors,  for  a  30mm  depth  of  burial 
simulation,  are  shown  in  Figures  16  through  20.  Figure  16  shows  that  there  was 
no  air  shock  wave  measured  by  PG1  for  the  30mm  depth  of  burial  simulation. 
Figure  17  shows  a  strange  wave  profile  in  the  shape  of  a  ramp  wave,  indicative 
of  a  supported  shock,  (which  should  not  be  the  case  here  as  the  lOkJ  is 
deposited  into  the  gold  hohlraum  within  a  few  nanoseconds),  for  the  air  pressure 
gauge  (PG2). 


Pressure  at  Pressure  Gauge  1  3cm  DOB 
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Figure  16.  Pressure  Wave  at  Air  Gauge  1  for  3cm  DOB 
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PrtttiM*  at  Ptataur*  Gauga  2  3cm  OOB 
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Figure  17.  Strange  Pressure  Wave  at  Air  Gauge  2  for  3cm  DOB 

Upon  initial  review  of  the  data,  it  was  unclear  how  the  pressure  profile  in 
Figure  17  translated  to  what  was  happening  mathematically  in  the  simulation. 
The  wave  profile  appears  to  indicate  a  pressure  wave  reflection  that  came  from 
the  boundary.  To  verify  that  this  was  indeed  the  case,  further  calculation  was 
conducted  in  which  an  additional  air  sensor  (PG3)  was  set  up  with  an  r-z 
coordinate  of  [25,140],  which  is  closer  to  the  boundary  of  the  setup,  to  see  if  the 
pressure  peak  time  of  arrival  was  shorter  than  that  in  PG2.  The  results,  shown  in 
Figure  18,  indicate  the  same  anomaly  in  the  pressure  wave  at  the  nominal  3rd  air 
pressure  gauge  (PG3).  The  peak  pressure  time  of  arrival  for  the  nominal  PG3 
sensor,  arriving  at  approximately  40  microseconds,  was  shorter  than  the  same 
anomalous  wave  in  Figure  17  for  PG2.  This  shows  that  the  data  is  consistent 
with  a  wave  that  reflected  from  the  boundary. 
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Pressure  at  Pressure  Gauge  3  3cm  DOB 


Figure  18.  Strange  Pressure  Wave  at  Nominal  Air  Pressure  Gauge  3 

In  order  to  eliminate  the  reflected  wave  noise  from  the  depth  of  burial 
calculations  the  set  up  of  the  problem  had  to  be  modified.  The  boundary 
conditions  of  the  problem  are  set  up  so  that  the  system  conserves  energy,  mass, 
etc.  This  issue  did  not  arise  in  the  height  of  burst  calculation  because  there  is  a 
single  material,  “air,”  at  the  boundary.  The  issue  arose  in  the  depth  of  burial 
calculations  because  the  boundary  has  two  different  materials,  both  Macor  and 
gas-mix.  To  correct  the  setup,  a  thin  layer  of  gas-mix  was  added  around  the 
Macor  object  in  the  input  file  and  the  simulations  were  run  again  (See  Appendix 
B).  Figure  19  shows  the  modified  depth  of  burial  setup.  When  performing  the 
pressure  analysis  on  the  modified  depth  of  burial  calculations,  the  noise  in  PG2 
vanished.  Additionally  the  peak  pressures  in  the  three  “ground”  sensors  in  the 
Macor  object  were  identical,  which  shows  that  the  modification  of  the  setup  did 
not  adversely  affect  the  peak  pressure  readings.  Figures  20  and  21  show  the 
pressure  wave  profiles  for  PG2  and  nominal  PG3,  which  show  that  little  to  no 
shock  was  transmitted  from  the  ground  to  the  gas-mix  (the  pressure  signal  is 
down  an  entire  order  of  magnitude  from  that  in  Figures  17  and  18),  as  expected. 


30 


Figure  19.  Modified  Depth  of  Burial  Set  Up 


x  10'°  Pressure  at  Pressure  Gauge  2  3cm  DOB 


Figure  20.  Pressure  at  PG2  with  Modified  Set  up  (no  noise  from  boundary) 
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Pressure  at  Pressure  Gauge  3  3cm  DOB 


Figure  21 .  Pressure  profile  for  nominal  pressure  gauge  3  (PG3) 


Pressure  at  Ground  Pressure  Gauge  1  3cm  DOB 


Figure  22.  Pressure  Wave  at  Ground  Gauge  1  for  3cm  DOB 
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Pressure  at  Ground  Pressure  Gauge  2  3cm  DOB 


Figure  23.  Pressure  Wave  at  Ground  Gauge  2  for  3cm  DOB 


Pressure  at  Ground  Pressure  Gauge  3  3cm  DOB 


Figure  24.  Pressure  Wave  at  Ground  Gauge  3  for  3cm  DOB 

As  one  would  intuitively  guess,  the  peak  pressures  calculated  in  the  Macor 
object  are  much  higher  for  the  depth  of  burial  shots.  Again  for  Figures  22  through 
24,  the  many  reflection  stress  waves  are  the  result  of  free  surface  reflections 
from  the  Macor/gas-mix  boundary. 
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VI.  ENERGY  COUPLING 


In  order  to  illustrate  energy  coupling  and  partitioning,  two  methods  were 
used.  The  first  method  demonstrates  the  coupling  by  looking  at  the  peak 
pressures  for  two  pressure  sensors  for  varied  heights  of  burst  and  depths  of 
burial.  Figure  25  shows  three  pressure  plots  for  three  different  heights  of  burst/ 
depth  of  burial  from  air  pressure  gauge  1  (PG1)  superimposed  on  one  another. 
Figure  26  shows  three  stress  plots  for  three  different  heights  of  burst/  depth  of 
burial  from  a  ground  sensor,  (GPG3),  superimposed  on  one  another.  These  two 
sensors,  PG1  and  GPG3,  were  chosen  because  both  are  located  at  a  distance  of 
10cm  from  the  point  of  detonation.  The  r-z  coordinates  for  PG1  are  (8.6cm,  5cm) 
and  the  r-z  coordinates  for  GPG3  are  (7.07cm,  7.07cm)  both  yielding  a  range  of 
10cm.  The  three  pressure  curves  correspond  to  two  heights  of  burst  and  one 
depth  of  burial.  The  plot  in  Figure  25  shows  that  the  air  blast  is  strongest  for 
detonations  that  occur  near  the  surface.  This  may  seem  counterintuitive.  The 
reason  that  the  surface  burst  yields  a  stronger  air  blast  is  that  the  ground  reflects 
much  of  the  pressure  wave,  and  the  shock  sensors  detect  the  combined 
pressure  of  both  the  initial  and  the  reflected  shock. 
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Figure  25.  Air  blast  is  strongest  for  surface  burst 
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As  reasonably  expected,  the  results  of  the  ground  pressure/stress  plots 
show  that  the  ground  shock  is  strongest  for  a  buried  burst. 


Figure  26.  Ground  shock  is  the  strongest  for  a  buried  burst 

The  second  method  used  to  illustrate  energy  coupling  and  partitioning  was 
to  calculate  the  equivalent  yield-coupling  factor.  This  was  done  by  plotting  the 
total  energy  in  the  Macor  block  as  a  function  of  time  for  each  height  of  burst  and 
depth  of  burial  calculation.  The  total  energy  transmitted  into  the  Macor  block  for 
the  height  of  burst  simulations  is  straightforward.  The  total  energy  in  the  Macor 
varies  inversely  with  the  height  of  burst.  At  the  theoretical  limit  of  an  infinite 
height  of  burst,  there  would  be  no  energy  transferred  into  the  Macor  block,  and 
the  only  energy  in  the  Macor  would  be  due  to  its  inherent  mass  energy.  Figure  27 
shows  the  energy  transmitted  into  the  Macor  block  as  a  function  of  time  for  each 
of  the  height  of  burst  calculations. 
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HOB  Energy  in  Macor 
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Figure  27.  Energy  in  Macor  block  for  several  HOB  shots 


As  expected,  the  energy  decreases  to  a  baseline  as  the  height  of  burst 
increases.  A  similar  relationship  should  also  hold  for  the  depth  of  burial  shots. 
At  the  theoretical  limit  of  an  infinitely  buried  depth  of  burst,  all  the  energy  will 
remain  within  the  Macor  and  none  transmitted  into  the  atmospheric  gas  mix. 
Figure  28  shows  the  energy  contained  within  the  Macor  block  as  a  function  of 
time  for  the  depth  of  burial  calculations. 
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Figure  28.  Energy  in  Macor  block  for  several  DOB  shots 
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The  total  energy  within  the  Macor  increases  from  5mm  depth  of  burst  to 
10mm  depth  of  burst,  yet  there  is  a  noticeable  difference  between  the  profiles  of 
the  5mm  and  10mm  DOB  energies  and  the  others.  As  the  depth  of  burial 
increases  beyond  20mm  DOB,  the  energy  begins  to  decline  with  depth  of  burst. 
At  first  glance,  it  may  seem  to  indicate  an  error  in  the  data,  as  the  energy  should 
increase  to  an  asymptote,  as  the  depth  of  burial  gets  deeper.  Upon  further 
investigation,  using  LLNL’s  visualization  application,  VisIT,  it  was  discovered  that 
for  both  the  5mm  and  10mm  depth  of  burial  runs  the  exploding  gas  escapes  the 
Macor  block  and  leaves  behind  a  crater.  The  point  at  which  the  Macor  gives  way 
corresponds  to  the  two  peaks  in  the  total  energy  curve  shown  above  for  both  the 
5mm  and  10mm  DOB  shots,  as  shown  in  Figure  28.  Figure  29  shows  a  split 
view  movie  frame,  with  material  on  the  right  and  pressure  on  the  left,  at  35 
microseconds  of  the  detonation.  The  breakout  of  the  gold  gas  from  the  Macor 
block  is  clear.  This  breakout  corresponds  to  the  peak  in  the  energy  versus  time 
plot  in  the  Macor  block  for  the  10mm  DOB  shot. 


Figure  29.  Split  view  of  Material  and  Pressure  at  35  microseconds  for  10mm 

DOB  shot 
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For  all  shots  below  10mm,  the  fireball  never  escapes  the  Macor  block,  but 
simply  bulges  the  “ground”  where  the  fireball  is  forming.  Figure  30  shows  that 
the  gold  gas  never  escapes  the  Macor  block  by  120  microseconds  for  a  20mm 
depth  of  burial  simulation. 


4.:  iJOOl 


Figure  30.  Split  view  of  Material  and  Pressure  at  120  microseconds  for  20mm 

DOB  shot 

As  demonstrated  in  Figures  30  and  31,  the  layer  of  Macor  surrounding  the 
gold  gas  becomes  thicker  as  the  depth  of  burial  increases. 
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Figure  31 .  Split  view  of  Material  and  Pressure  at  120  microseconds  for  30mm 

DOB  shot 

Thus,  for  20mm  and  lower  depths  of  burial,  the  energy  contained  within 
the  gold  gas  fireball  was  added  to  the  total  energy  of  the  Macor  block,  and  with 
this  addition,  the  total  energy  asymptotes  as  expected.  The  baseline  energy  of 
the  Macor,  due  to  its  mass  energy,  was  subtracted  from  each  to  show  that  the 
energy  asymptotes  to  the  lOkJ  value  of  total  deposited  energy.  Figure  32 
illustrates  the  results  of  this  corrected  calculation  of  energy  in  the  Macor  block. 
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DOB  Energy  in  Macor 
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Figure  32.  Corrected  Plot  of  Energy  in  Macor  block  for  several  DOB  shots 


With  the  corrected  energy-time  plots  for  both  height  of  burst  and  depth  of 
burial  runs,  the  peak  energy  values  were  plotted  as  a  function  of  depth  of  burial 
in  Figure  33  in  order  to  construct  the  more  familiar  energy  coupling  curve. 


Energy  Coupling  Cuive  (Lab  Frame) 


Figure  33.  Lab  Frame  Energy  Coupling  Curve 
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In  this  plot,  the  ordinate  corresponds  to  the  energy  coupling  yield  factor 
(ECYF),  which  is  simply  the  peak  energy  in  the  Macor  divided  by  the  total  input 
energy  of  10kJ.  A  zero  ECYF  would  correspond  to  an  infinite  high  height  of  burst 
detonation  and  an  ECYF  of  one  would  correspond  to  an  infinitely  deep  depth  of 
burial  detonation.  The  abscissa  indicates  depth  of  burial  where  negative  values 
along  this  axis  correspond  to  the  heights  of  burst  for  above-ground  detonations. 
One  key  aspect  of  these  data  is  that  they  are  scalable.  In  order  to  scale  up  from 
the  lab  frame,  in  units  of  millimeters  and  kilojoules,  the  cube  root  scaling  law, 
explained  earlier  in  Chapter  II  sub-section  B,  “Scaling,”  was  used  to  convert  the 
straight  depth  of  burial  to  a  scaled  depth  of  burial  with  units  of  distance  per  yield 
energy^^^.  The  conversion  factor  shown  below  takes  advantage  of  the  cube  root 
scaling  law  noted  by  Glasstone  and  Dolan. 

Conversion  factor: 

\kJ  =  2.383*10  (megaton) 

(lAJ)*''  =(2.383*10  =61.99*10 '’MY*'" 

\mm  =  0.001m 

1mm  /  lAJ*'"  =  0.00 Im  /  (6 1 .99  *  1 0  )  =  1 6. 1 3m  /  MY*'' 

Each  of  the  depth  of  burial/  height  of  burst  was  scaled  up  to  a  scaled 
depth  of  burial/  scaled  height  of  burst  using  the  conversion  factor.  How  each 
distance  was  scaled  is  shown  below,  using  the  5mm  depth  of  burial  case  as  an 
example. 

simulation  depth  of  burial/  yield'^^  *  conversion  factor  =  scaled  depth  of  burial 

5mm  / 1  * (1 6. 1 3m  / )  /  (1mm  /  =  31 36m  I Mt'^ 

Table  5  shows  the  results  of  similar  calculations  for  each  depth  of  burial 
ranging  from  5mm  to  50  mm  with  a  10kJ  yield. 
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mm  /10kJ''^ 

m  /MT"'^ 

5 

37.36 

10 

74.7 

20 

149.45 

30 

224.18 

40 

298.9 

50 

373.64 

Table  5.  Scaled  Depth  of  Burial 


Figure  34  shows  the  energy-coupling  curve  scaled  to  the  world  frame. 


Energy  Coupling  Curve  Scaled  Up  to  Real  World  Frame 


Figure  34.  World  Frame  Energy  Coupling  Curve 

The  data  can  now  be  mapped  to  earlier  calculations  of  explosions  in 
various  hard  rock  media  made  by  several  researchers  and  summarized  by  L. 
Glenn  at  LLNL  in  Figure  35.  Those  calculations  are  focused  in  the  regime  where 
there  is  strong  coupling  between  the  ground  and  air  shock. 
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Energy  Coupled  into  the  Ground  is  a  Strong  Function 
of  Depth  of  Burial  (DOB) 
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Figure  35.  Energy  Coupling  Curve  for  Hard  Rock  From  L. Glenn,  2002,  p.16 


First,  the  upper  limit  of  the  data  fit  was  calculated  by  using  the  equation 
from  the  depth  of  burial  data  fit  in  Figure  35,  setting  the  equivalent  yield  coupling 
factor  (ECYF)  to  unity,  and  calculating  the  scaled  depth  of  burial  that 
corresponds  to  full  coupling,  i.e.,  no  energy  transmitted  into  air  shock. 

1  =  0.22  *  scaled  depth  of  burial  (SDOB)°  ^® 
ln(l  /  0.22)  =  0.36  *  \n{SDOB) 

\n{SDOB)  =  (1  /  0.36)  *  ln(l  /  0.22) 

SDOB  =  ^  67,08[m 

This  shows,  according  to  the  data  fit,  that  all  shots  conducted  at  a  scaled 
depth  of  burial  greater  than  67.08  m/  represent  a  situation  in  which  all  of 
the  energy  would  be  fully  coupled  into  the  ground.  Using  the  scaled  depth  of 
burial  points  from  Table  2  and  the  data  fit  formulas  summarized  by  L.  Glenn,  the 

ECYF  data  points  from  the  numerical  simulation  was  plotted  against  the  coupling 
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curve  of  the  hard  rock  data  in  order  to  show  the  correlation  with  the  Macor 
simulation  model.  Figure  36  shows  that  within  the  limits  of  the  data  fit  model, 
both  the  simulation  and  the  hard  rock  model  data  fit  agree  within  a  margin  of  5%. 
The  red  triangular  points  are  from  the  simulation  data  for  the  various  HOB/DOB 
calculations.  The  green  curve  is  a  plot  of  the  two  formulas  from  the  hard  rock 
data  equations. 


Matched  Energy  Coupling  Curve 


Figure  36.  Matched  Coupling  Curve 
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Table  6  shows  the  EYCF  for  each  of  the  12  simulation  calculations.  The 
EYCF  was  calculated  by  taking  the  peak  energies,  from  Figures  27  and  32,  for 
the  corresponding  height  of  burst  and  depth  of  burial  calculations  respectively, 
and  then  dividing  that  value  by  the  lOkJ  input  energy.  So,  for  example,  in  the 
10mm  DOB  calculation  the  maximum  energy  contained  within  the  Macor  was 
9.4kJ  which  corresponds  to  the  peak  of  the  10mm  curve  in  Figure  32.  This  value 
of  9.4kJ  was  then  divided  by  the  lOkJ  yield  to  result  in  an  ECYF  of  0.94. 


Depth  of  Burial/  Height  of  Burst 

Equivalent  Yield  Coupling  Factor 

5mm 

.85 

10mm 

.94 

20mm 

.99 

30mm 

.99 

40mm 

.99 

50mm 

.99 

-5mm  (Height  of  Burst) 

.012 

-10mm 

.006 

-20mm 

.003 

-30mm 

.0025 

-40mm 

.002 

-50mm 

.0018 

Table  6.  Equivalent  Yield  Coupling  Factor  from  simulation  data 


Table  7  shows  a  similar  data  chart,  the  difference  being  that  the  scaled 
depths  of  burial/  heights  of  burst  (shown  in  Table  5)  were  used  along  with  the  two 
equations  from  the  hard  rock  data  for  both  the  height  of  burst  and  depth  of  burial 
to  calculate  the  EYCF. 
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Scaled  DOB/  Scaled  HOB  [m/MT"'^] 

Equivalent  Yield  Coupling  Factor 

37.36 

.81 

74.7 

SDOB  larger  than  67.08  ^1 

149.45 

SDOB  larger  than  67.08  ^1 

224.18 

SDOB  larger  than  67.08  ^1 

298.9 

SDOB  larger  than  67.08  ^1 

373.64 

SDOB  larger  than  67.08  ^1 

-37.36  (Height  of  Burst) 

.012  (**outside  data  range  yet  agrees) 

-74.7 

.011  (outside  data  range) 

-149.45 

.010  (outside  data  range) 

-224.18 

.010  (outside  data  range) 

-298.9 

.010  (outside  data  range) 

-373.64 

.009  (outside  data  range) 

Table  7.  Equivalent  Yield  Coupling  Factor  from  live  data  fit  curves 


The  range  in  which  both  the  simulation  data  and  the  hard  rock  data 
overlap  is  small,  yet  within  that  range  both  the  simulation  and  hard  rock  data 
agree  to  within  5%.  Even  slightly  beyond  the  limits  of  the  curve  fit  from  the  hard 
rock  data  both  the  simulation  and  it  agree,  as  noted  in  the  5mm  HOB  (-37.36 
m/MT^'^  scaled  HOB),  where  both  the  simulation  and  the  hard  rock  data  fit  predict 
the  exact  same  ECYF  of  .012. 
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VII.  CONCLUSION 


Scaled  simulation  and  computer  modeling  of  nuclear  weapon  effects  will 
play  a  significant  role  in  the  maintenance  of  the  SBSSP.  The  National  Ignition 
Facility  at  LLNL  provides  a  unique  opportunity  to  investigate  nuclear  weapon 
effects  and  the  physics  of  nuclear  weapons  in  support  of  the  achievement  of 
goals  put  forth  by  the  NNSA  and  the  SBSSP.  This  work  provides  information  in 
evaluating  the  hydrodynamic  response  of  materials  to  nuclear  detonations  at  or 
near  the  earth’s  surface,  which  is  a  major  component  of  the  SPSSP  agenda. 
The  work  presented  in  this  thesis  serves  two  purposes:  first,  to  assist  in  the 
planning  of  the  placement  of  shock  sensors  within  a  Macor  object  that  will  be 
used  as  a  platform  for  upcoming  NIF  experiments  designed  to  study  the  energy 
coupling/partitioning  of  a  low  yield  nuclear  detonation  at  or  near  the  earth’s 
surface;  and  second,  to  develop  a  basic  set  of  computational  data  for 
development  of  a  coupling  curve.  The  coupling  curve  derived  here,  though  partly 
based  on  an  idealistic  buried  device,  has  proven  to  be  very  consistent  with 
calculations  of  previous  explosions  in  various  hard  rock  media  supported  by 
limited  test  data.  The  coupling  curve  derived  here  and  the  supporting 
calculations  serve  as  a  basic  set  of  data  for  further  investigation  of  buried  nuclear 
devices.  Using  the  work  done  in  this  thesis,  opportunities  for  further  research 
include  adding  the  effect  of  radiation  to  the  simulations  and  altering  the  design  of 
the  Macor  object  to  more  closely  model  an  inhomogeneous  material. 
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APPENDIX  A 


Appendix  A  is  the  input  file  for  a  10mm  height  of  burst  calculation.  Many 
of  the  original  marker  files  have  been  removed  for  the  purpose  of  saving  paper. 
This  should  help  clarify  particular  portions  of  the  text,  specifically  the  topics  of 
boundary  conditions  and  equations  of  state.  This  input  file  was  called  into  the 
GEODYN  executable  and  the  results  analyzed  with  both  MATLAB  and  VisIT. 


%  units  are  specified  in  millimeters,  milligrams,  microseconds  (GPa,  km/s) 

%%%%%%%%%%%%%%%%%%%%%%% 

%%%%  DEFINITIONS  %%%%%% 
%%%%%%%%%%%%%%%%%%%%%%% 

HOB=10  %  height  of  burst 

%%%%%%%%%%%%%%%%%%%%%%% 

%%  GLOBAL  PARAMETERS  %% 
%%%%%%%%%%%%%%%%%%%%%%% 
max_step  =  1 00000 
stop_time  =  200 


%%%%%%%%%%%%%%%%%%%%%%% 
%%%  DOMAIN  GEOMETRY  %%% 
%%%%%%%%%%%%%%%%%%%%%%% 
geometry:( 

%  cylindrical  (r-z)  coordinates 
coord_sys  =  1 

%  low  boundaries  of  (r,z)  in  domain 
probjo  =  [0  -120] 

%  high  boundaries  of  (r,z)  in  domain 
prob_hi  =  [  160  160] 

) 


%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%  MESH/OUTPUT  %%%%% 
%%%%%%%%%%%%%%%%%%%%%%% 
amr:( 

%  number  of  elements  in  (r,z)  on  base  grid 
n_cell  =  [80  140] 
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%  verbose  output 
V  =  1 

%  verbose  timestep  output 
v_step  =  1 

%  maximum  level  of  refinement 
maxjevel  =  2 

%  regrid  interval  (timesteps) 
regridjnt  =  2 
%  maximum  grid  size 
max_grid_size  =  32 
%  refinement  ratio  (per  level) 
ref_ratio  =  [44444] 

%  number  of  extra  cells  around  a  grid  to  refine  (per  level) 
n_error_buf  =  [12  1  1  1] 

%  root  name  of  restart  files 
check_file  =  sesAucheck 
%  restart  file  interval  (steps) 
checkjnt  =100 

%  root  name  of  plot  files 
plot_file  =  sesAuPlot 

%  plot  file  interval  (steps)  [cannot  be  specified  with  plot_per] 

%plotJnt  =100 

%  plot  file  period  (microseconds)  [cannot  be  specified  with  plotjnt] 
plot_per  =  1 

%  plot  variables  (see  RULESETS) 

plot_vars  =  [  pressure  density  color  vr  vz  sound_speed  internal_energy 
bulking_porosity  pressure_cutoff  friction_slope  Til  T22  T33  T12  temperature  ] 

%  marker  specifications  [  name  (l)agrange/(e)ulerian  r  z  ] 


Markers  =  [ 

[ pgaugel 

e  86.0 

50.0] 

[ pgauge2 

e  130.0 

75.0] 

[ pgauge3 

e  140.0 

25.0] 

[ gpgauge2 

1  0.5 

-70.0  ] 

[ gpgaugel 

1  35.0 

-35.0  ] 

[ gpgauge3 

1  70.0 

-70.0  ] 

] 

%  marker  output  interval  (steps) 
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markerjnt  =  1 

%  marker  file  write  interval  (steps) 
marker_dumpjnt  =  100 
%  marker  variables  (see  RULESETS) 

marker_vars  =  [  m_pressure  m_density  m_vel_mag  m_vr  m_vz  m_T11  m_T22 
m_T33  m_T12  m_color  m_temperature  m_internal_energy  m_soundspeed 
m_plastic_strain  m_Tmin  m_rdisp  m_zdisp  ] 

) 


%%%%%%%%%%%%%%%%%%%%%%% 

%%%  PHYSICS/SOLVER  %%%% 

%%%%%%%%%%%%%%%%%%%%%%% 

hyp:( 

%  number  of  materials 
n_materials=5 
%  verbosity 

V  =  0 

%  minimum  timestep  below  which  to  halt 
dt_cutoff  =  0.00001 
%  Courant  stability  factor 
cfl  =  .3 

%  initial  scale  factor  for  timestep 
init_shrink  =  .01 

%  maximum  timestep  increase  per  cycle 
change_max  =1.1 
%  gravitational  acceleration 
gravity  =  0. 

%  interval  to  sum  integrated  quantities 
sumjnterval  =100 

%  integrated  quantity  variables  (see  RULESETS) 

integrated_quantities  =  [  sum_mass_source  sum_energy_source 

sum_mass_source2  sum_energy_source2  sum_mass_source0 

sum_energy_source0  ] 

%  write  integrated  quantities  file  on  restart 
%  i  n  teg  rate_on_resta  rt=  1 

) 


%%%%%%%%%%%%%%%%%%%%%%% 
%  PROBLEM  CONDITIONS  %% 
%%%%%%%%%%%%%%%%%%%%%%% 
Prob:( 

%  materials  (see  MATERIALS) 
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materials=[  gas_mix_LEOS  Au_LEOS  Macor  Air_LEOS  Fe_LEOS  ] 

%  boundary  condition  types 

%%  0  =  Interior/Periodic  2  =  Outflow  %% 

%%  1  =  Inflow  3  =  Symmetry  %% 

bc:lo:type  =  [31] 
bc:hi:type  =  [11] 

%  inflow  boundary  specifications 
%%  material  is  0-indexed,  density  in  g/cc,  energy  in  kJ/g 
bc:lo;1:( 
material=3 
rho=0.0001 
e=0.52 

) 

bc:hi;0:( 

material=3 

rho=0.0001 

e=0.52 

) 

bc:hi:1:( 

material=3 

rho=0.0001 

e=0.52 

) 


%  Specifies  that  displacement  variables  should  be  stored 
celldata<CellData>:( 
ops:[ 

<CellDisplacement>=0 

] 

) 


%%%%%%%%%%%%%%%%%%%%%%% 

%%%%  CONFIGURATION  %%%% 
%%%%%%%%%%%%%%%%%%%%%%% 

Config:( 

%  List  of  region  names 

%%  Specified  order  that  regions  are  painted  in. 

%%  Last  region  fills  remaining  space. 

Regions  =  [  source  air_cone  steel_cone  ceramic 
ring  atmosphere  ceramic_shell  air  ] 

%  Region  specifications 

%%  material  is  0-indexed,  density  in  g/cc,  energy  in  kJ/g 
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Region:( 

%  source 

%%  10kJ  in  2mg  of  material,  sphere  of  diameter  2mm 
source:( 

%  expanded  iron 
material  =  1 
rho  =  0.4774648 
e  =  5e+3 
file  =  circle512 

%  no  scaling,  translate  in  +y  by  height  of  burst 
Transform  =  [  scale  1.1.  translate  0.  '$:HOB'  ] 

) 

%  air  inside  cone 
air_cone:( 

%  air 

material  =  3 
rho  =  0.0001 
e  =  0.52 

%  r-z  pairs  specifying  polygon  (r-z  projection  of  cone) 
Points  =  [ 

[0  '$;HOB+1000'] 

['0.6+1000/6'  '$:HOB+1001'] 

[0.6  '$:HOB+1'  ] 

[0.  '$:HOB'  ] 

] 


%  steel  cone 
steel_cone;( 

%  iron  at  reference  density  and  energy 
material  =  4 

%  r-z  pairs  specifying  polygon  (r-z  projection  of  outer  cone) 
Points  =  [ 

['0.4+1000/6'  '$:HOB+1001'] 

[0.4  '$:HOB+1'  ] 

[0.6  '$;HOB+1'  ] 

[1.4  '$:HOB+1.8'  ] 

['1.4+1000/6'  '$:HOB+1001.8'] 

] 


%  macor 
ceramic:( 

%  macor  at  reference  density 
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material  =  2 
e  =  0.0000405 
Points  =  [ 

[0  -100] 
[100  -100] 
[100  0.06] 
[0  0.06] 

] 


%  macor  ring 
ring:( 

%  macor  at  reference  density 
material  =  2 
e  =  0.0000405 


Points  =  [ 

[  100.5 

-10] 

[  149.5 

-10] 

[  150. 

0.06] 

[  100. 

0.06] 

] 


%  gas-mixture  atmosphere 
atmosphere:( 

%  gas  mixture 
material  =  0 
rho=0.00153 
temperature  =  310 
%  composition  of  mixture 
"qq[0]"  =  0.65 
"qq[ir  =  0.20 

"qq[2]"  =  0.10 
"qq[3]"  =  0.05 
Points  =  [ 

[0  -100] 

[150  -100] 

[150  150] 

[0  150] 

] 

) 

%  macor  shell 
ceramic_shell:( 

%  macor  at  reference  density 
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material  =  2 


e  =  0.0000405 
Points  =  [ 


[0 

-105] 

[  155 

-105] 

[  155 

155] 

[0 

155] 

[0 

150] 

[  150 

150] 

[  150 

-100] 

[0 

-100] 

) 

%  remaining  air 
air:( 

%  air 

material  =  3 
rho  =  0.0001 
e  =  0.52 

) 

) 

) 


%%%%%%%%%%%%%%%%%%%%%%% 

%%%  REFINEMENT  RULES  %% 
%%%%%%%%%%%%%%%%%%%%%%% 

%  Specifies  refinement  condition 
refinement<Tagger>:( 

%  refine  on  material  interfaces  and  pressure  gradients 
logic<SymbolFunction>:( 
f="(interfaces|pressure)" 
variable="interfaces, pressure" 

) 

%  refinement  variables  (see  RULESETS) 
derive=[ref_color  ref_p  ] 

%  material  interface  detection 
%%  refine  when  "color"  is  not  an  integer 
interfaces<FabFunction>:( 
comp=[ref_color] 

f=abs_remainder_offset  parameters=[1  .e-8] 

) 
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%  pressure  gradient  detection 

%%  refine  when  relative  difference  (pA  -  pB)/(min(pA,pB,pMin))  >  0.3 
%%  where  pMin  =  0.0002GPa 
pressure<FabFunction>:( 
comp=[ref_p] 

f=relative_gradient_offset  parameters=[0.0002  .3] 

) 


) 


%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%  RULESETS  %%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%% 


%  plot  variables  (rules  must  match  names) 
pressure<DeriveRuleSet>:( 
models=[  StressTensor  Energies  ] 
rules:[ 

<Rule>:(var=p  averager=bulk_modulus_average) 
<Rule>:(var=rho  averager=volume_average) 
<Rule>:(var=material_number  averager=volume_average) 
<Rule>:(var="v[0]" ) 

<Rule>:(var="v[1]" ) 

<Rule>:(var=c) 

<Rule>:(var=e) 

<Rule>:(var=bulking_porosity  ) 
<Rule>:(var=pressure_cutoff ) 

<Rule>:(var=friction_slope) 

<Rule>:(var="T[T11]" ) 

<Rule>:(var="T[T22]" ) 

<Rule>:(var="T[T33]" ) 

<Rule>:(var="T[T12]" ) 

<Rule>:(var=temperature) 

] 

) 

%  marker  variables  (rules  must  match  names) 
m_pressure<DeriveRuleSet>:( 
models=[StressTensor  Energies  ] 
rules;[ 

<Rule>:(var=p  averager=bulk_modulus_average) 
<Rule>:(var=rho  averager=volume_average) 
<Rule>:(var=v_mag) 

<Rule>:(var="v[0]" ) 

<Rule>:(var="v[1]" ) 
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<Rule>:(var="T[T11]") 

<Rule>:(var="T[T22]") 

<Rule>:(var="T[T33]") 

<Rule>:(var="T[T12]") 

<Rule>:(var=material_number  averager=volume_average) 
<Rule>:(var=temperature) 

<Rule>:(var=e) 

<Rule>:(var=c) 

<Rule>:(var=plastic_strain) 

<Rule>:(var=Tminhistory) 

<Rule>:(var="displacement[0]") 

<Rule>:(var="displacement[1]") 

] 


%  integrated  quantities  variables  (rules  must  match  names) 
sum_mass_source<DeriveRuleSet>:( 
models  =  [  Energies  ] 
rules;[ 

<Rule>:(var=rho  averager=volume_sum  materials=[1]) 
<Rule>:(var=e_tot  averager=mass_sum  materials=[1]) 
<Rule>:(var=rho  averager=volume_sum  materials=[2]) 
<Rule>:(var=e_tot  averager=mass_sum  materials=[2]) 
<Rule>:(var=rho  averager=volume_sum  materials=[0]) 
<Rule>:(var=e_tot  averager=mass_sum  materials=[0]) 

] 

) 

%  refinement  variables 
ref_color<DeriveRuleSet>:( 
models=[EOS] 
rules:[ 

<Rule>:(var=material_number  averager=volume_average) 
<Rule>:(var=p  averager=bulk_modulus_average) 

] 

) 


%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%  MATERIALS  %%%%%% 
%%%%%%%%%%%%%%%%%%%%%%% 
gas_mix_LEOS<Material>:  ( 
models:  [ 

<NoPorosityVariable>  =  0 
<RhoSolid>  =  0 
<N_GammaAveraged_EOS>:( 
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EOS;[ 

<LEOS>  :(  use_energy_offset=1  material_number=100  )  %  qq[0]  -  Ne 
<LEOS>  :(  use_energy_offset=1  material_number=180  )  %  qq[1]  -  Ar 
<LEOS>  :(  use_energy_offset=1  material_number=360  )  %  qq[2]  -  Kr 
<LEOS>  :(  use_energy_offset=1  material_number=540  )  %  qq[3]  -  Xe 

] 

) 

<N_LimitFraction>  =  0 
<SphericalStressTensor>  =  0 
<ZeroShearModulus>  =  0 
<Smooth_PressureCutOff>:  (pmin=0) 

<ZeroFrictionSlope>  =  0 
<Tminhistory>  =  0 

] 

) 

Air_LEOS<Material>:  ( 
models:  [ 

<NoPorosityVariable>  =  0 
<RhoSolid>  =  0 

<LEOS>  :(  material_number=2260  ) 

<NoUpdate>  =  0 
<SphericalStressTensor>  =  0 
<ZeroShearModulus>  =  0 
<Smooth_PressureCutOff>:  (pmin=0) 

<ZeroFrictionSlope>  =  0 
<Tminhistory>=0 

] 

) 

Fe_LEOS<Material>:  ( 
models:  [ 

<NoPorosityVariable>  =  0 
<RhoSolid>  =  0 

<LEOS>:  (  material_number=260  ) 

<NoUpdate>  =  0 
<SphericalStressTensor>  =  0 
<ZeroShearModulus>  =  0 
<Smooth_PressureCutOff>:  (pmin=0) 

<ZeroFrictionSlope>  =  0 
<Tminhistory>=0 

] 

) 

Au_LEOS<Material>:  ( 
models:  [ 
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<NoPorosityVariable>  =  0 
<RhoSolid>  =  0 

<LEOS>:  (  material_number=2700  leos_name=/usr/gapps/data/eos/sesame  ) 

<NoUpdate>  =  0 

<SphericalStressTensor>  =  0 

<ZeroShearModulus>  =  0 

<Smooth_PressureCutOff>:  (pmin=0) 

<ZeroFrictionSlope>  =  0 
<Tminhistory>=0 

] 

) 


Macor<Material>:  ( 
models:  [ 

<ElasticPorosityVariable>  =  0 
<RhoSolid>  =  0 
<MieGruneisen_EOS>:  ( 
rhoO  =  2.52 
cO  =4.58 
s1  =1.4 
s2  =0. 
s3  =0. 

g  =  1- 

b  =  0. 

) 

<HistoryDependent_PressureCutOff>:  (  pmin=-0.008  ) 
<PoreElasticity>:  (  a=0.  b=1.  ) 

<PseudoCapYieldStrength1>:( 
scale_factor=1 . 

Y_c=0.345 

Y_t=0.1 

tau_dam=0. 00001 
P_hard_rate=5. 

P_hard_exp=0.5 

Yrat=0.98 

cap0=0.5 

cap00=0.5 

cap_power=2. 

strain_hard_eps=0.02 

strain_tofail=0.0001 

soft_rate=5. 

rate_function<SymbolicFunction>:  ( function="(x+1  .)''0.12" ) 
residual_function<SymbolicFunction>:  ( function="1  .*x" ) 

) 
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<PseudoCapCompaction>:( 

poro0=0.001 

mu_c0=10. 

slope0=0.1 

bulking_k1<SymbolicFunction>:  ( function="0.05*x" ) 
tau_comp=0 

) 

<RateDependentFunctionPlasticUpdatewithCorrection>:( 
tau_of_p<Const>:  ( f=0  ) 

) 

<ElasticUpdate1>  =  0 

<YieldLodeAngleTerm>:  (  q20  =  0.5  function<Const>:(f=1) ) 

<Stress>  =  0 

<ConstantPoissonRatio>:  (  nu=0.29  ) 

<PseudoCapBulkingwithCorrection>  =  0 

<PseudoCapBulkingCorrection>:  ( 
k1=0.5 
k=0.5 

) 

<FrictionSlope>:  ( y0=1.5 ) 

] 
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APPENDIX  B 


Appendix  B  is  the  input  file  for  the  modified  30mm  depth  of  burial 
calculation.  Many  of  the  original  marker  files  have  been  removed  for  the  purpose 
of  saving  paper.  The  original  goal  of  the  modified  DOB  set  up  was  to  have  a 
layer  of  air  around  both  the  Macor  and  the  gas-mix,  but  without  a  solid  boundary 
between  the  two  gases  the  gas-mix  expanded  to  fill  the  problem  space.  This 
setup  worked  also  as  it  still  left  a  single  material  at  the  problem  boundary,  which 
eliminated  the  reflected  wave  noise,  as  shown  in  Figures  17  and  18.  This  input 
file  was  called  into  the  GEODYN  executable  and  the  results  analyzed  with  both 
MATLAB  and  VisIT. 

%  units  are  specified  in  millimeters,  milligrams,  microseconds  (GPa,  km/s) 

%%%%%%%%%%%%%%%%%%%%%%% 

%%%%  DEFINITIONS  %%%%%% 
%%%%%%%%%%%%%%%%%%%%%%% 

FIOB=-30  %  depth  of  burst 

%%%%%%%%%%%%%%%%%%%%%%% 

%%  GLOBAL  PARAMETERS  %% 

%%%%%%%%%%%%%%%%%%%%%%% 
max_step  =  1 00000 
stop_time  =  200 


%%%%%%%%%%%%%%%%%%%%%%% 
%%%  DOMAIN  GEOMETRY  %%% 
%%%%%%%%%%%%%%%%%%%%%%% 
geometry:( 

%  cylindrical  (r-z)  coordinates 
coord_sys  =  1 

%  low  boundaries  of  (r,z)  in  domain 
probjo  =  [0  -120] 

%  high  boundaries  of  (r,z)  in  domain 
prob_hi  =  [  160  160] 

) 


%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%  MESH/OUTPUT  %%%%% 
%%%%%%%%%%%%%%%%%%%%%%% 
amr:( 
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%  number  of  elements  in  (r,z)  on  base  grid 
n_cell  =  [80  140] 

%  verbose  output 
V  =  1 

%  verbose  timestep  output 
v_step  =  1 

%  maximum  level  of  refinement 
maxjevel  =  2 

%  regrid  interval  (timesteps) 
regridjnt  =  2 
%  maximum  grid  size 
max_grid_size  =  32 
%  refinement  ratio  (per  level) 
ref_ratio  =  [44444] 

%  number  of  extra  cells  around  a  grid  to  refine  (per  level) 
n_error_buf  =  [12  1  1  1] 

%  root  name  of  restart  files 
check_file  =  chk 
%  restart  file  interval  (steps) 
checkjnt  =  200 

%  root  name  of  plot  files 
plot_file  =  pit 

%  plot  file  interval  (steps)  [cannot  be  specified  with  plot_per] 

%plotJnt  =100 

%  plot  file  period  (microseconds)  [cannot  be  specified  with  plotjnt] 
plot_per  =  1 

%  plot  variables  (see  RULESETS) 

plot_vars  =  [  pressure  density  color  vr  vz  sound_speed  internal_energy 
bulking_porosity  pressure_cutoff  friction_slope  Til  T22  T33  T12  temperature  ] 

%  marker  specifications  [  name  (l)agrange/(e)ulerian  r  z  ] 


Markers  =  [ 

[  pgaugel 

e  86.0 

50.0] 

[  pgauge2 

e  130.0 

75.0] 

[  pgauge3 

e  25.0 

140.0] 

[  gpgauge2 

1  0.5 

-70.0  ] 

[  gpgaugel 

1  35.0 

-35.0  ] 

[  gpgauge3 

1  70.0 

-70.0  ] 
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] 

%  marker  output  interval  (steps) 
markerjnt  =  1 

%  marker  file  write  interval  (steps) 
marker_dumpjnt  =  100 
%  marker  variables  (see  RULESETS) 

marker_vars  =  [  m_pressure  m_density  m_vel_mag  m_vr  m_vz  m_T11  m_T22 
m_T33  m_T12  m_color  m_temperature  m_internal_energy  m_soundspeed 
m_plastic_strain  m_Tmin  m_rdisp  m_zdisp  ] 

) 


%%%%%%%%%%%%%%%%%%%%%%% 

%%%  PHYSICS/SOLVER  %%%% 

%%%%%%%%%%%%%%%%%%%%%%% 

hyp:( 

%  number  of  materials 
n_materials=4 
%  verbosity 

V  =  0 

%  minimum  timestep  below  which  to  halt 
dt_cutoff  =  0.00001 
%  Courant  stability  factor 
cfl  =  .2 

%  initial  scale  factor  for  timestep 
init_shrink  =  .01 

%  maximum  timestep  increase  per  cycle 
change_max  =1.1 
%  gravitational  acceleration 
gravity  =  0. 

%  interval  to  sum  integrated  quantities 
sumjnterval  =100 

%  integrated  quantity  variables  (see  RULESETS) 

integrated_quantities  =  [  sum_mass_source  sum_energy_source 

sum_mass_source2  sum_energy_source2  sum_mass_source0 

sum_energy_source0  ] 

%  write  integrated  quantities  file  on  restart 
%  i  n  teg  rate_on_resta  rt=  1 

) 


%%%%%%%%%%%%%%%%%%%%%%% 
%  PROBLEM  CONDITIONS  %% 
%%%%%%%%%%%%%%%%%%%%%%% 
Prob:( 
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%  materials  (see  MATERIALS) 

materials=[  gas_mix_LEOS  Au_LEOS  Macor  Air_LEOS  ] 

%  boundary  condition  types 

%%  0  =  Interior/Periodic  2  =  Outflow  %% 

%%  1  =  Inflow  3  =  Symmetry  %% 

bc:lo:type  =  [31] 
bc:hi:type  =  [11] 

%  inflow  boundary  specifications 
%%  material  is  0-indexed,  density  in  g/cc,  energy  in  kJ/g 
bc:lo:1:( 
material=3 
rho=0.0001 
e=0.52 

) 

bc:hi:0:( 

material=3 

rho=0.0001 

e=0.52 

) 

bc:hi:1:( 

material=3 

rho=0.0001 

e=0.52 

) 


%  Specifies  that  displacement  variables  should  be  stored 
celldata<CellData>:( 
ops:[ 

<CellDisplacement>=0 

] 

) 


%%%%%%%%%%%%%%%%%%%%%%% 
%%%%  CONFIGURATION  %%%% 
%%%%%%%%%%%%%%%%%%%%%%% 
Config:( 

%  List  of  region  names 

%%  Specified  order  that  regions  are  painted  in. 
%%  Last  region  fills  remaining  space. 

Regions  =  [  source  ceramic  atmosphere  ] 

%  Region  specifications 
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%%  material  is  0-indexed,  density  in  g/cc,  energy  in  kJ/g 
Region:( 

%  source 

%%  10kJ  in  2mg  of  material,  sphere  of  diameter  2mm 
source:( 

%  expanded  gold 
material  =  1 
rho  =  0.4774648 
e  =  5e+3 
file  =  circle512 

%  no  scaling,  translate  in  +y  by  height  of  burst 
Transform  =  [  scale  1.1.  translate  0.  '$:HOB'  ] 

) 

%  macor 
ceramic:( 

%  macor  at  reference  density 
material  =  2 
e  =  0.0000405 
Points  =  [ 

[0  -110] 

[150  -110] 

[150  0.06] 

[0  0.06] 

] 


%  gas-mixture  atmosphere 
atmosphere:( 

%  gas  mixture 
material  =  0 
rho=0.00153 
temperature  =  310 
%  composition  of  mixture 
"qq[0]"  =  0.65 
"qq[ir  =  0.20 

"qq[2]"  =  0.10 
"qq[3]"  =  0.05 
Points  =  [ 

[0  0.06] 

[0  150] 

[150  150] 

[150  0.06] 

] 

) 
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%  air  boundary 
air:( 

%  air  boundary 
material  =  3 
e  =  0.52 
rho=0.0001 

) 


) 


) 


%%%%%%%%%%%%%%%%%%%%%%% 

%%%  REFINEMENT  RULES  %% 
%%%%%%%%%%%%%%%%%%%%%%% 

%  Specifies  refinement  condition 
refinement<Tagger>:( 

%  refine  on  material  interfaces  and  pressure  gradients 
logic<SymbolFunction>:( 
f="(interfaces|pressure)" 
variable="interfaces, pressure" 

) 

%  refinement  variables  (see  RULESETS) 
derive=[ref_color  ref_p  ] 

%  material  interface  detection 
%%  refine  when  "color"  is  not  an  integer 
interfaces<FabFunction>:( 
comp=[ref_color] 

f=abs_remainder_offset  parameters=[1  .e-8] 

) 

%  pressure  gradient  detection 

%%  refine  when  relative  difference  (pA  -  pB)/(min(pA,pB,pMin))  >  0.3 
%%  where  pMin  =  0.0002GPa 
pressure<FabFunction>:( 
comp=[ref_p] 

f=relative_gradient_offset  parameters=[0.0002  .3] 

) 


) 


%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%  RULESETS  %%%%%%% 
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%%%%%%%%%%%%%%%%%%%%%%% 


%  plot  variables  (rules  must  match  names) 
pressure<DeriveRuleSet>:( 
models=[  StressTensor  Energies  ] 
rules:[ 

<Rule>:(var=p  averager=bulk_modulus_average) 
<Rule>:(var=rho  averager=volume_average) 
<Rule>:(var=material_number  averager=volume_average) 
<Rule>:(var="v[0]" ) 

<Rule>:(var="v[1]" ) 

<Rule>:(var=c) 

<Rule>:(var=e) 

<Rule>:(var=bulking_porosity  ) 
<Rule>:(var=pressure_cutoff ) 

<Rule>:(var=friction_slope) 

<Rule>:(var="T[T11]" ) 

<Rule>:(var="T[T22]" ) 

<Rule>:(var="T[T33]" ) 

<Rule>:(var="T[T12]" ) 

<Rule>:(var=temperature) 

] 

) 

%  marker  variables  (rules  must  match  names) 
m_pressure<DeriveRuleSet>:( 
models=[StressTensor  Energies  ] 
rules:[ 

<Rule>:(var=p  averager=bulk_modulus_average) 
<Rule>:(var=rho  averager=volume_average) 
<Rule>:(var=v_mag) 

<Rule>:(var="v[0]" ) 

<Rule>:(var="v[1]" ) 

<Rule>:(var="T[T11]") 

<Rule>:(var="T[T22]") 

<Rule>:(var="T[T33]") 

<Rule>:(var="T[T12]") 

<Rule>:(var=material_number  averager=volume_average) 
<Rule>:(var=temperature) 

<Rule>:(var=e) 

<Rule>:(var=c) 

<Rule>:(var=plastic_strain) 

<Rule>:(var=Tminhistory) 

<Rule>:(var="displacement[0]") 

<Rule>:(var="displacement[1]") 
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] 

) 

%  integrated  quantities  variables  (rules  must  match  names) 
sum_mass_source<DeriveRuleSet>:( 
models  =  [  Energies  ] 
rules:[ 

<Rule>:(var=rho  averager=volume_sum  materials=[1]) 
<Rule>:(var=e_tot  averager=mass_sum  materials=[1]) 
<Rule>:(var=rho  averager=volume_sum  materials=[2]) 
<Rule>:(var=e_tot  averager=mass_sum  materials=[2]) 
<Rule>:(var=rho  averager=volume_sum  materials=[0]) 
<Rule>:(var=e_tot  averager=mass_sum  materials=[0]) 


] 

) 

%  refinement  variables 
ref_color<DeriveRuleSet>:( 
models=[EOS] 
rules:[ 

<Rule>:(var=material_number  averager=volume_average) 
<Rule>:(var=p  averager=bulk_modulus_average) 

] 

) 


%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%  MATERIALS  %%%%%% 
%%%%%%%%%%%%%%%%%%%%%%% 
gas_mix_LEOS<Material>:  ( 
models:  [ 

<NoPorosityVariable>  =  0 
<RhoSolid>  =  0 
<N_GammaAveraged_EOS>:( 

EOS:[ 

<LEOS>  :(  use_energy_offset=1  material_number=100  )  %  qq[0]  -  Ne 
<LEOS>  :(  use_energy_offset=1  material_number=180  )  %  qq[1]  -  Ar 
<LEOS>  :(  use_energy_offset=1  material_number=360  )  %  qq[2]  -  Kr 
<LEOS>  :(  use_energy_offset=1  material_number=540  )  %  qq[3]  -  Xe 

] 

) 

<N_LimitFraction>  =  0 
<SphericalStressTensor>  =  0 
<ZeroShearModulus>  =  0 
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<Smooth_PressureCutOff>:  (pmin=0) 
<ZeroFrictionSlope>  =  0 
<Tminhistory>  =  0 

] 


Au_LEOS<Material>:  ( 
models:  [ 

<NoPorosityVariable>  =  0 
<RhoSolid>  =  0 

<LEOS>:  (  material_number=2700  leos_name=/usr/gapps/data/eos/sesame  ) 

<NoUpdate>  =  0 

<SphericalStressTensor>  =  0 

<ZeroShearModulus>  =  0 

<Smooth_PressureCutOff>:  (pmin=0) 

<ZeroFrictionSlope>  =  0 
<Tminhistory>=0 

] 

) 

Air_LEOS<Material>:( 
models:  [ 

<NoPorosityVariable>=0 

<RhoSolid>=0 

<LEOS>:(material_number=2260) 

<NoUpdate>=0 
<SphericalStressTensor>=0 
<ZeroShearModulus>=0 
<Smooth_PressureCutOff>:  (pmin=0) 

<ZeroFrictionSlope>=0 

<Tminhistory>=0 

] 

) 


Macor<Material>:  ( 
models:  [ 

<ElasticPorosityVariable>  =  0 
<RhoSolid>  =  0 
<MieGruneisen_EOS>:  ( 
rhoO  =  2.52 
cO  =4.58 
s1  =1.4 
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s2  =0. 
s3  =0. 

g  =  i- 

b  =  0. 

) 

<HistoryDependent_PressureCutOff>:  (  pmin=-0.008  ) 
<PoreElasticity>:  (  a=0.  b=1.  ) 

<PseudoCapYieldStrength1>:( 
scale_factor=1 . 

Y_c=0.345 

Y_t=0.1 

tau_dam=0. 00001 
P_hard_rate=5. 

P_hard_exp=0.5 

Yrat=0.98 

cap0=0.5 

cap00=0.5 

cap_power=2. 

strain_hard_eps=0.02 

strain_tofail=0.0001 

soft_rate=5. 

rate_function<SymbolicFunction>:  ( function="(x+1  .)''0.12" ) 
residual_function<SymbolicFunction>:  ( function="1  .*x" ) 

) 

<PseudoCapCompaction>:( 

poro0=0.001 

mu_c0=10. 

slope0=0.1 

bulking_k1<SymbolicFunction>:  ( function="0.05*x" ) 
tau_comp=0 

) 

<RateDependentFunctionPlasticUpdatewithCorrection>:( 
tau_of_p<Const>:  ( f=0  ) 

) 

<ElasticUpdate1>  =  0 

<YieldLodeAngleTerm>:  (  q20  =  0.5  function<Const>:(f=1) ) 
<Stress>  =  0 

<ConstantPoissonRatio>:  (  nu=0.29  ) 
<PseudoCapBulkingwithCorrection>  =  0 
<PseudoCapBulkingCorrection>:  ( 
k1=0.5 
k=0.5 

) 

<FrictionSlope>:  (  y0=1 .5  ) 

] 

) 
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